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Abstract. - A new multiple-relaxation-time lattice Boltzmann scheme for compressible flows with 
arbitrary speciflc heat ratio and Prandtl number is presented. In the new scheme, which is based 
on a two-dimensional 16-discrete- velocity model, the kinetic moment space and the corresponding 
transformation matrix are constructed according to the seven-moment relations associated with 
the local equilibrium distribution function. In the continuum limit, the model recovers the com- 
pressible Navier-Stokes equations with flexible specific-heat ratio and Prandtl number. Numerical 
experiments show that compressible flows with strong shocks can be simulated by the present 
model up to Mach numbers Ma ~ 5. 



Introduction. — Over the last two decades, the Lat- 
tice Boltzmann (LB) method has become a prominent tool 
in computational fluid dynamics [IHS] , with a broad spec- 
trum of complex flows applications, ranging from mul- 
tiphase flows [IMj, magnetohydrodynamics [THTU]. flows 
through porous media [llH13j and thermal fluid dynam- 
ics [14l[T5], to cite but a few. Being based on a low- 
Mach expansion of Maxwell-Boltzmann local equilibria, 
standard LB models are typically used for the simula- 
tion of quasi-incompressible flows. Given the great impor- 
tance of shock-wave dynamics in many fields of physics 
and engineering [T6l|T7] , constructing LB models for high- 
speed compressible flows with shocks has been attempted 
since the early days of LB research. However, current 
LB versions for compressible flows are still subject to at 
least one of the following constraints: low Mach num- 
ber, fixed Prandtl number and fixed specific heat ratio. 
Broadly speaking, to date, there are two LB approaches 
for the simulation of compressible flows. The first is based 
on the Single-Relaxation-Time (SRT) Bhatnagar-Gross- 
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Krook (BGK) approximation [T^, where the local equilib- 
rium distribution function f'^'^ is calculated from a trun- 
cated Taylor expansion of the flow velocity, where i is the 
index of discrete velocity. Models along this line can be 
roughly categorized into two groups, the standard LB and 
the Finite Difference (FD) LB. Examples in the first group 
are referred to the works of Alexander and Chen et al. [19] , 
Sun, et al. [10], Yan, et al. [H], and so on. In the FDLB 
formulation, space and time discretizations are indepen- 
dent, and to overcome the limit of low Mach number, one 
has to confine the LB to be a solver of Euler or Navier- 
Stokes (NS) equations. In such a case, the discrete /f' is 
best seen as a local attractor of the coUisional relaxation 
process, leading to the desired form of the macroscopic mo- 
ments. Examples are referred to the works of Tsutahara, 
et al [HI13], Xu, et al [H], Shu, et al [IS] and He, et al 
[55]. The application of BGK also leads to a fixed Prandtl 
number. In most existing models, the specific heat ratio 
is fixed to a nonrealistic constant because only the trans- 
lational degrees of freedom are taken into account in the 
internal energy. The second way consists of resorting to 
the original scattering-matrix formulation of the LB equa- 
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tion |27j . whose optimized version is nowadays known as 
Multiple- Relaxation-Time (MRT) LB scheme [25]. In the 
MRT-LB formulation, the collision step is first calculated 
in the kinetic moment space (KMS), i.e. the space spanned 
by the kinetic moments defined by the projection of the 
distribution function onto a basis set of N Hcrmite polyno- 
mials. Subsequently, the streaming step is performed back 
in the discrete velocity space spanned by a suitable set of 
N discrete velocities v^. In contrast to the SRT model, the 
MRT version caters for more adjustable parameters and 
degrees of freedom. The relaxation rates of the various ki- 
netic moments due to particle collisions may be adjusted 
independently. This overcomes some obvious deficiencies 
of the SRT model, such as a fixed Prandtl number. The 
low-Mach number constraint is generally caused by the 
numerical instability problem. To this regard, extensive 
efforts have been made in the past years, such as the en- 
tropic method [MS], the fix-up scheme Flux- 
limiters f33[ and dissipation |241l34j techniques. In many 
cases with low-speed fiows, it has been shown [35] that the 
MRT LB model also offers enhanced numerical stability. 
As far as we know, nearly all existing MRT models have 
been focussing on the standard LB for isothermal systems 
without strong shocks. 

Recently, we presented a thermal MRT FDLB model 
for high-speed compressible flows [55]. The MRT model 
uses 16 discrete velocities, as suggested by Kataoka and 
Tsutahara [23] in the SRT version. In that MRT model, 
the KMS and corresponding transformation matrix are 
constructed according to the group representation theory; 
equilibria of the nonconserved moments are chosen so as to 
recover the compressible Navier-Stokes equations through 
the Chapman-Enskog analysis. Neither the formulation, 
nor the simulation procedure are by any means related to 
a truncated Taylor-expansion of local equilibrium distribu- 
tion function f^"^. Consequently, the discrete local equi- 
librium f^'^ remains Maxwellian regardless of the value of 
the Mach number. However, the model is restricted to a 
nonphysical value of the specific-heat-ratio, 2. For conve- 
nience of description, we will refer to that model as model 
/. 

In this letter, we formulate a LB model II for compress- 
ible flows with arbitrary speciflc-heat-ratio. At variance 
with model I, at high-Mach number, the local equilibria 
of model II depart from a Maxwellian. Nevertheless, the 
method is able to simulate both low and high Mach num- 
bers regimes (up to Ma ^ 5). 

Brief review of the MRT LB modeL According 
to the main strategy of the MRT-LB scheme, the differen- 
tial form of the MRT-LB equation read as follows: 

^+--^ = -M,7^S,,(A-/f), (1) 

where is the discrete particle velocity, i = 1,. . . ,N, N 
is the number of discrete velocities, the subscript a indi- 
cates X or y. The variable t is time, Xa is the spatial coor- 
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Fig. 1: Distribution of Vi for the discrete velocity model. 

dinate. The matrix S = MSM~^ = diag{si, S2, • • • , sat) 
is the diagonal relaxation matrix, fi and fi are the par- 
ticle distribution function in the velocity space and the 
kinetic moment space respectively, fi = niijfj, rriij is 
an element of the matrix M. Obviously, the mapping 
between KMS and velocity space is defined by the lin- 
ear transformation M, i.e., f = Mf, f = M~^f, where 
the bold-face symbols denote N-dinicnsional column vec- 
tors, e.g., f = (/i,/2,--- JnV, f = (/i,/2,--- JnV, 
M = {nil, 1112, ■ ■ ■ jTOat)'^, mi = {mii,mi2, ■ ■ ■ ,m^N)- 
is the equilibrium value of the moment /,;. The equation 
reduces to the usual lattice BGK equation if all the relax- 
ation parameters are set to be a single relaxation time r, 
namely S = -l-I, where I is the identity matrix. 

Construction of energy- conserving MRT LB 

model. — Wc use the two-dimensional discrete veloc- 
ity model by Kataoka and Tsutahara (sec Fig. 1), which 
reads as follows: 




eye: (±1,0), for 1 < I < 4, 

eye: (±6,0), for 5 < i < 8, 

y2(±l,±l), for 9 < i < 12, 

^(±1,±1), forl3<^<16, 



where eye indicates the cyclic permutation. The 16 ve- 
locities are grouped into four energy levels, u^, where 
V = 1,2,3,6. In this model, besides the translational de- 
grees of freedom, a parameter rji is introduced, in order to 
describe the (&— 2) extra-degrees of freedom corresponding 
to molecular rotation and/or vibration, where rji = 5/2 for 
z = 1, • • • ,4, and iii = for 1 = 5, • • • , 16. Details are re- 
ferred to the original publication [23] . In our simulations, 
the time evolution is based on the usual first-order upwind 
scheme, while space discretization is performed through a 
Lax-Wcndroff scheme. 

The moment representation provides a convenient way 
to express coUisional relaxation on different time scales for 
the various kinetic moments. However, the choice of the 
kinetic moments is not unique. They can be identified 
with physical quantities, such as density, momentum, en- 
ergy, momentum and energy fluxes, and so on. Inspired 
by a close tie with tensor Hermite polynomials, many re- 
searchers choose the moments according to the monomials 
of Cartesian components of the discrete velocities v™v^ , 
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m, n = 0, 1, 2, • • • [3S]. In this Letter, the KMS and equi- 
libria of the moments are chosen according to the seven- 
moment relations (Equations (5a-5g) in Reference [131), 
associated with the local equilibrium distribution function 
/j'^'' (Equation (8) in Reference [13]), which is a polynomial 
of the flow velocity up to the third order. The incorpora- 
tion of the parameter b permits to tune the specific-heat- 
ratio, -f = {b + 2)/b. 

The right-hand-sides of the seven equations indicate 
seven monomials, 1, Via, {vf^ + ijf), ViaVip, {vfp + vD^ia, 
ViaVipVi-^, {vf^ + rif)viaVip. In order to make the choice 
of moments more transparent, we construct the transfor- 
mation matrix through a simple combination of the seven 
monomials above (sec Appendix for details). As a result, 
at least eight of these moments carry an explicit physical 
meaning: /i = mufi is the fluid density; /2 — m2ifi and 
/a = m^ifi correspond to the x— and y~ components of 
momentum (mass flux); /4 = mufi corresponds to the to- 
tal energy mode; /g = rriQifi and f-^ = m-afi correspond 
to the diagonal and off-diagonal components of the stress 
tensor; /g — mgifi and /g = rngifi correspond to the x— 
and y— components of the total energy flux. 

At variance with previous MRT models [351 [37], the 
transformation matrix M used in this work is not based 
upon a Gram-Schmidt orthogonalization procedure. This 
is because such a procedure does not lead to a conserved 
energy, as required for a consistent thermal and compress- 
ible model. The moments can be divided into two groups. 
The flrst group consists of those which are locally con- 
served in the collision process, i.e. fi = f^''. There are 
four such conserved moments, density p, momenta jx, jy, 
and energy e in two dimensions, which are denoted by 
/i, /2, /s, /4 mentioned above, respectively. The second 
group consists of the non-conserved ones, i.e. fi ^ f^''. 
The equilibria of the non-conserved moments are function- 
als of the conserved ones. For the collision process, one 
has maximum freedom in the construction of the equilib- 
rium functions of the non-conserved moments, provided 
the basic principles of physics are satisfied (conservations 
of mass, momentum and energy). In this model, the cor- 
responding equilibrium distribution functions in the KMS 
f^"^ are chosen according to the seven-moment relations, 
too (see Appendix). By using the Chapman-Enskog ex- 
pansion [33] on the two sides of the LB equation, the final 
NS equations for compressible fluids read as follows: 



dp , djx djy 



djx 
dt 



dt dx 
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dP 
dx 



(2a) 



d_ 

dy 



Jxjy 




dy ^ bs5 



+ (2c) 

se dx S7 



| + |:[(^ + 2P)Wp] + |[(e + 2P)Vp] 

^d ,pRT,,b ,,^dT ,dux 2^, 

„ d ,pRT,,b ^dT 



(2d) 



dx 0S5 sg dy sj 



(2b) 



where P = pRT, e = bpRT + j'^/ P is twice the total 
energy, s^y = (dxUy + dyU^), D = {d^u^ + dyUy), and 
d = (dxUx — dyUy), S5, sg, S7 are related to viscosity, sg, sg 
are related to heat conductivity. The relaxation param- 
eters are not completely independent, since the isotropy 
constraints impose coupling between some of these relax- 
ation parameters, e.g., sg = sg. Since the elements of 
S represent the inverse of the relaxation time for f to 
its equilibrium f^"* in KMS, and the values of p, jx^ jy, 
and e are conserved in the relaxation process, the val- 
ues of si,S2, 53,54 can be fine-tuned independently. The 
remaining relaxation parameters must belong to the in- 
terval [0, l/dt] for reasons of linear stability. In the limit 
S5 = 56 = 57 = 5g = 59 = 5, the above NS equations 
reduce to the SRT compressible NS equations. 

Numerical Simulation and Analysis. In this 
section we study the following problems using the new 
MRT LB model: measure of the sound speed, one- and 
two-dimensional Riemann problems, Richtmyer-Meshkov 
instability. We work in a frame where the constant R = 1. 

(a) Speed of sound 

In order to validate the proposed model, the sound 
speed is calculated and compared to the theoretical value 
Cs = y/lT [H]. A one-dimensional tube is used. In 
the tube, a plate divides the fluids, that share the same 
temperature and feature a small difference in pressure 
{Pi/P2 = 1.0 + 10^^°) at the initial stage. The pressure 
ratio is so small that expansion and compression waves 
propagate at the sound speed in both directions when the 
plate is removed. The position of the pressure jump is 
measured to calculate the speed of sound. The results at 
various temperature levels are shown in Fig. 2, and the 
calculated values are found to agree well with the theoret- 
ical ones. 

(h) One- dimensional Riemann problem 

Here, the two-dimensional model is used to solve the 
one-dimensional Riemann problem, so as to visually show 
the advantages of this model as compared to its SRT ver- 
sion. Consider the Riemann problem with initial condition 



p-3 



Chen, Xu, Zhang, Li and Succi 



Theoretical value y=1 -667 

Theoretical value y=1 .400 

o FDLBM simulation 



1.0 

temperature T 



Fig. 2: Calculation results and theoretical values of sound 
speed. 



given by 

{p,ui,u2,T)\l = (5.99924,19.5975,0.0,76.8254), 
(p, ui,u2,T)\r = (5.99242, -6.19633, 0.0, 7.69222), 

where subscript "L" and "R" indicate the left and right 
macroscopic variables of discontinuity. The numerical 
and exact solutions for two different specific-heat ratios 
7 = 7/5 and 7 = 5/3 at time t — 0.13 are shown in Fig. 3, 
where the common parameters for SRT and MRT simula- 
tions are dx ^ dy ^ 0.002, dt = 10^^. The relaxation time 
in SRT is T = 10~^, while the collision parameters in MRT 
are S5 = 5x 10^, sg = lO"', and other values of s are lO'^. In 
the X direction, fi = f^"^ is set on the boundary nodes be- 
fore the disturbance reaches the two ends. The equilibrium 
distribution function f^"^ in velocity space can be obtained 
via the linear transformation, = M~^f'^'?. Along the y 
direction, periodic boundary conditions are adopted. The 
red triangle and red circle symbols correspond to MRT 
simulation results with 7 = 7/5 and 7 ~ 5/3, respectively; 
the green triangle and green circle symbols correspond to 
SRT simulation results with 7 = 7/5 and 7 = 5/3, respec- 
tively; the solid lines represent the exact solutions for the 
case 7 = 7/5, and the dashed lines represent the exact 
solutions for 7 = 5/3. We find that the oscillations at 
the discontinuity are weaker in MRT simulations than in 
their SRT countrepart. This shows that the problematic 
"wall-heating" phenomenon is weaker in the MRT model 
than in its SRT version. 

(c) Two-dimensional Riemann problem 

Compared with the relatively simple 1-D configurations, 
the 2-D Riemann problem consists of a plethora of geomet- 
ric wave patterns that pose a challenge for the simulation. 
In two-dimensions, the Riemann problem consists of four 
uniform states in each quadrant, (p, ui, W2, T)(x, y, 0) — 
[pi, uii, U2i, Ti), i = 1, 2, 3, 4, where i denotes the ith quad- 
rant. We solve the 2-D Riemann problem with initial data 
as follows: 

pi = 1.5,wii = 0,it2i = 0,ri = 1; 
P2 = 0.5323,ui2 = 1.206,U22 = 0,r2 = 0.3/0.5323; 
P3 = 0.138,ui3 = 1.206,U23 = 1.206,r3 = 0.029/0.138; 
P4 = 0.5323,ui4 = 0,U24 = 1.206,r4 = 0.3/0.5323. 
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Fig. 3: The numerical and exact solutions for the one- 
dimensional Riemann problem at time t = 0.13. 

Figure 4 shows the density contours of the 2-D Riemann 
problem by MRT LB simulation at time t — 0.12, where 
the parameters are dx = dy = 0.002, dt = 10~^, sg = sq = 
2 X lO'', other values of s are 3 x 10'^. The Mach number of 
the second and fourth quadrants are 1.358, and the third 
quadrant is 3.1444. While the SRT version is found to 
fail, the MRT LB successfully recovers the main features 
observed in earlier computations [5^140) . 

(d) Richtmyer - Meshkov instability 

Studies of the shock-induced Richtmyer - Meshkov 
(RM) instability are of great importance from both view- 
points of fundamental research and practical applications. 
For example, the RM instability occurs in Incrtial Con- 
fined Fusion, Scramjet engine, supersonic and hypersonic 
combustion, and others. In addition, it also plays an im- 
portant role in some natural phenomena, such as Super- 
nova explosions, formation of salt domains, and others. 
Although most practical problems are three-dimensional, 
2D studies are also of recognized value, as they offer useful 
insights into the basic physics of these problems. The sim- 
ulation of RM instability still remains a challenging topic 
to this day. 

The initial configuration of our simulation is illustrated 
in Fig. 5(a), where an interface with sinusoidal per- 
turbation separates two different fluids and a incident 
shock wave, with the Mach number 2.5, traveling from 
the right side, hits the interface immediately. The com- 
putational domain is a rectangle with length 0.36 and 
height 0.1, which is divided into 360 x 100 mesh-cells. 
The initial sinusoidal perturbation at the interface is: 
X = 0.24 4- 0.005 X sin(7r/2 + iO-Ky), where the cycle of 
initial perturbation is 0.05, and the amplitude is 0.005. 
The initial conditions arc as follows: 

{p,ui,u2,p)i = (0.1358,0,0,1), 

{p,Ul,U2,p)rn = (1,0,0, 1), 

{p,ui,U2,p)r = (3.33333,-2.07063,0,7.125), 
where the subscripts I, m, r indicate the left, middle, right 
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Fig. 4: (Color online) Density contour lines of the 2-D Riemann 
problem in MRT LB simulation at time t = 0.12. From purple 
to gray, the density value increases. 
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Fig. 5: (Color online) Snapshots of RM instability: density 
contour lines at times t = 0, t = 0.02, t = 0.074, t = 0.104, 
respectively. From blue to red the color corresponds to the 
increasing of density. 



regions of the whole domain. The following boundary con- 
ditions are imposed: (1) inflow at the right side; (2) re- 
flecting condition at the left boundary, and (3) periodic 
boundary conditions are applied at the top and bottom 
boundaries. The reflecting boundary condition implies 
that the x component of the fluid velocity on the bound- 
ary is the reverse of the one of the mirror-image interior 
point. The parameters are as follows: 7 = 1.4, = 10~^, 
.S5 = 2 X 10^, and 10^ for the others. 

When the shock wave passes the interface from the right, 
a reflected rarefaction wave to the right and a transmis- 
sion wave to the left, are generated, and the perturbation 
amplitude decreases with the interface motion to the left 
(Figure 5 (b)). Then, the peak and valley of initial inter- 
face invert, the heavy and light fluids gradually penetrate 
into each other as time goes on, the light fluid "rises" to 
form a bubble and the heavy fluid "falls" to generate a 
spike ( shown in Figure 5 (c)). The transmission wave 
reaches the solid wall on the left and reflects to the right, 
encounters the interface again, converts into a transmis- 
sion wave penetrating into the heavy fluid zone and a re- 
flection wave back to the light fluid. Before the second 
shock wave reaches the interface, the interface continues to 
move, and the disturbance develops at a lower growth rate. 
When the second shock reaches the interface, the instabil- 
ity growth rate of the disturbance signiflcantly enhances. 
The disturbance of the interface continues to grow, even- 
tually forming a mushroom shape (Figure 5 (d)). The sim- 
ulation results show a satisfactory agreement with those 
of other numerical simulations j41U42j and experiments 
[321111] ■ In the experimental study of Puranik et al. [13] 
the shock wave travels from light (air) to heavy (CO2) 
medium, while the reverse is true in our simulation. As a 
result, interface reversals are observed in our simulation. 
The same mechanism is also mentioned in Ref. j43) (see 
Fig. 11 in [131), which is a typical characteristic of an 
incident shock from the heavier to the lighter material. 



Even though the region with light medium is a cylindrical 
bubble in their experiments, while it is a rectangle with 
distortions in our simulation, the experiments on the inter- 
action of a shock wave with a single cylindrical bubble [Hj 
also show high similarity with our simulation results, both 
in terms of the reversal mechanism and of the mushroom 
shape. 

Conclusions and remarks. — In this Letter, we have 
proposed a MRT-FDLB model for compressible flows with 
flexible specific heat ratio and Prandtl number, as appro- 
priate for most applications. Different transport coeffi- 
cients, such as viscosity and heat conductivity, are related 
to different collision parameters, thereby allowing a sep- 
arate control of the different transport processes. The 
new scheme has been validated for a series of one and 
two-dimensional numerical benchmarks, always showing 
satisfactory agreement with theoretical results and previ- 
ous numerical work. Three-dimensional extensions of the 
present scheme appear conceptually straightforward and 
shall make the object of future work. 
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Nos. 2009A0102005, 2009B0101012], National Natural 
Science Foundation of China [under Grant Nos. 10775018, 
10702010]. F. Chen and Y. Li acknowledge support of Na- 
tional Basic Research Program of China [under Grant No. 
2007CB815105]. 

Appendix. Construction of the transformation 
matrix and f^"^ . — In order to construct the transforma- 
tion matrix in KMS, we take the monomial ViaVip as an ex- 
ample. Three possibilities arise: (a) a = /3 = x, ViaVip = 
vfx, (b) a = ^ = y, WjaWj/? = vfy, (c) a = x, (3 = y, 
ViaVip = VixViy. "(a)-l-(b)" gives {vf^+vfy), "(a)-(b)" gives 
(^L ~ ^iy)- this way, the transformation matrix can 
be composed as follows: M = (mi,m2, • • • ,toi6)"^, where 
mii = 1, TO2i = m^i = Viy, = vf^ + 
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VixViy, mst = Vi^{vf^ + vfy +rif), mgi = Viy{vf^ + 

v% + ), miQ, = Vi^{vf^ + vfy); mil, = Viy{vf^ + 

vfy), mi2, = Vi^{vf^ - vfy), mi3, = Viy{vf^ - 

vfy), mu, = {vf^ + vfy){vf^ + vfy + Tjf), TOisi = 

Viy {vf^ +Vfy+7]f), mi6, = [vf^ - vfy ) ivf^+ vfy + 7]f), 

where i = 1, ■ ■ ■ , 16. 

The corresponding equilibrium distribution functions in 
KMS are as follows: f^'' = p, f^'^ = j^, /g' = 

01 - jDIp. fj' - J.Jy/P, fl' = (e + 
2P)jWp, /r = (e + 2P)jy/P. /To = (4^ + fJP + 
fy/p)j./p, m - (4P + jIIp + Jllp)jylp, Al = 
(2P + pjp - 3l/p)3./p, /Tl = (-2P + 31/ p 
3llp)3ylp, III = 2(5 + 2)pi?2T2 + (6 + h)RT[3l + 
3l)lp + 01 + jI?Ip\ /Tl = [(& + 4)P + 01 + 
J1)/P]J.JWP^ /Tl = [(^'+4)P+01+J,')/P]01-J,^)/P'- 
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